Cplex
Padova, 27/03/2024
Cari studenti,
vi riporto qui sotto alcuni frammenti di codice C (non ancora funzionanti: sono da
testare/estendere e debuggare!) sviluppati "alla lavagna" come esempio di codifica
in C degli algoritmi studiati durante il mio corso on-line di RO2 a.a. 2019-2020.
Buon divertimento!
MF
#include <cplex.h>
...
typedef struct {
//input data
int nnodes;
double *xcoord;
double *ycoord;
int integer_costs; // = 1 for integer costs (rounded distances), 0 otherwise
// parameters
int model_type;
double timelimit; // overall time limit, in sec.s
char input_file[1000]; // input file
//global data
double zbest; // value of the best sol. available
} instance; //----JUST DONE---
...
void print_error(const char *err) {
printf("\n\n ERROR: %s \n\n", err);
fflush(NULL);
exit(1);
}
double dist(int i, int j, instance *inst){
double dx = inst->xcoord[i] - inst->xcoord[j];
double dy = inst->ycoord[i] - inst->ycoord[j];
if ( !inst->integer_costs ) return sqrt(dx*dx+dy*dy);
int dis = sqrt(dx*dx+dy*dy) + 0.499999999; // nearest integer
return dis+0.0;
}
< prototipi di altre funzioni >
...
/**************************************************************************************************************************/
int TSPopt(instance *inst){
// open CPLEX model
int error;
CPXENVptr env = CPXopenCPLEX(&error);
CPXLPptr lp = CPXcreateprob(env, &error, NULL);
build_model(inst, env, lp);
// Cplex's parameter setting
// ...
if ( CPXmipopt(env,lp) ) print_error("CPXmipopt() error");
// use the optimal solution found by CPLEX
int ncols = CPXgetnumcols(env, lp);
double *xstar = (double *) calloc(ncols, sizeof(double));
//CPXgetx get the solution x* of our instance
if ( CPXgetx(env, lp, xstar, 0, ncols-1) ) print_error("CPXgetx() error");
for ( int i = 0; i < inst->nnodes; i++ )
{
for ( int j = i+1; j < inst->nnodes; j++ )
{
if ( xstar[xpos(i,j,inst)] > 0.5 ) printf(" ... x(%3d,%3d) = 1\n", i+1,j+1);
}
}
free(xstar);
// free and close cplex model
CPXfreeprob(env, &lp);
CPXcloseCPLEX(&env);
return 0; // or an appropriate nonzero error code
}
int xpos(int i, int j, instance *inst){
if ( i == j ) print_error(" i == j in xpos" );
if ( i > j ) return xpos(j,i,inst);
int pos = i * inst->nnodes + j - (( i + 1 ) * ( i + 2 )) / 2;
return pos;
}
void build_model(instance *inst, CPXENVptr env, CPXLPptr lp){
int izero = 0;
char binary = 'B';
char **cname = (char **) calloc(1, sizeof(char *)); // (char **) required by cplex...
cname[0] = (char *) calloc(100, sizeof(char));
// add binary var.s x(i,j) for i < j
for ( int i = 0; i < inst->nnodes; i++ )
{
for ( int j = i+1; j < inst->nnodes; j++ )
{
sprintf(cname[0], "x(%d,%d)", i+1,j+1); // x(1,2), x(1,3) ....
double obj = dist(i,j,inst); // cost == distance
double lb = 0.0;
double ub = 1.0;
if ( CPXnewcols(env, lp, 1, &obj, &lb, &ub, &binary, cname) ) print_error(" wrong CPXnewcols on x var.s");
if ( CPXgetnumcols(env,lp)-1 != xpos(i,j, inst) ) print_error(" wrong position for x var.s");
}
}
// add degree constr.s
int *index = (int *) malloc(inst->nnodes * sizeof(int));
double *value = (double *) malloc(inst->nnodes * sizeof(double));
// add the degree constraints
for ( int h = 0; h < inst->nnodes; h++ ) // degree constraints
{
double rhs = 2.0;
char sense = 'E'; // 'E' for equality constraint
sprintf(cname[0], "degree(%d)", h+1);
int nnz = 0;
for ( int i = 0; i < inst->nnodes; i++ )
{
if ( i == h ) continue;
index[nnz] = xpos(i,h, inst);
value[nnz] = 1.0;
nnz++
}
if ( CPXaddrows(env, lp, 0, 1, nnz, &rhs, &sense, &izero, index, value, NULL, &cname[0]) ) print_error(" wrong CPXaddrows [degree]");
}
free(value);
free(index);
if ( VERBOSE >= -100 ) CPXwriteprob(env, lp, "model.lp", NULL);
free(cname[0]);
free(cname);
}
#define DEBUG // da commentare se non si vuole il debugging
#define EPS 1e-5
/*********************************************************************************************************************************/
void build_sol(const double *xstar, instance *inst, int *succ, int *comp, int *ncomp) // build succ() and comp() wrt xstar()...
/*********************************************************************************************************************************/
{
#ifdef DEBUG
int *degree = (int *) calloc(inst->nnodes, sizeof(int));
for ( int i = 0; i < inst->nnodes; i++ )
{
for ( int j = i+1; j < inst->nnodes; j++ )
{
int k = xpos(i,j,inst);
if ( fabs(xstar[k]) > EPS && fabs(xstar[k]-1.0)) > EPS ) print_error(" wrong xstar in build_sol()");
if ( xstar[k] > 0.5 )
{
++degree[i];
++degree[j];
}
}
}
for ( int i = 0; i < inst->nnodes; i++ )
{
if ( degree[i] != 2 ) print_error("wrong degree in build_sol()");
}
free(degree);
#endif
*ncomp = 0;
for ( int i = 0; i < inst->nnodes; i++ )
{
succ[i] = -1;
comp[i] = -1;
}
for ( int start = 0; start < inst->nnodes; start++ )
{
if ( comp[start] >= 0 ) continue; // node "start" was already visited, just skip it
// a new component is found
(*ncomp)++;
int i = start;
int done = 0;
while ( !done ) // go and visit the current component
{
comp[i] = *ncomp;
done = 1;
for ( int j = 0; j < inst->nnodes; j++ )
{
if ( i != j && xstar[xpos(i,j,inst)] > 0.5 && comp[j] == -1 ) // the edge [i,j] is selected in xstar and j was not visited before
{
succ[i] = j;
i = j;
done = 0;
break;
}
}
}
succ[i] = start; // last arc to close the cycle
// go to the next component...
}
}
**** LAZY CONTRAINTS ****
Ex: MZT formulation with directed-arc variables x_ij and x_ji --> xpos_compact(i,j,inst)
...
int izero = 0;
int index[3];
double value[3];
// add lazy constraints 1.0 * u_i - 1.0 * u_j + M * x_ij <= M - 1, for each arc (i,j) not touching node 0
double big_M = inst->nnodes - 1.0;
double rhs = big_M -1.0;
char sense = 'L';
int nnz = 3;
for ( int i = 1; i < inst->nnodes; i++ ) // excluding node 0
{
for ( int j = 1; j < inst->nnodes; j++ ) // excluding node 0
{
if ( i == j ) continue;
sprintf(cname[0], "u-consistency for arc (%d,%d)", i+1, j+1);
index[0] = upos(i,inst);
value[0] = 1.0;
index[1] = upos(j,inst);
value[1] = -1.0;
index[2] = xpos_compact(i,j,inst);
value[2] = big_M;
if ( CPXaddlazyconstraints(env, lp, 1, nnz, &rhs, &sense, &izero, index, value, cname) ) print_error("wrong CPXlazyconstraints() for u-consistency");
}
}
// add lazy constraints 1.0 * x_ij + 1.0 * x_ji <= 1, for each arc (i,j) with i < j
rhs = 1.0;
char sense = 'L';
nnz = 2;
for ( int i = 0; i < inst->nnodes; i++ )
{
for ( int j = i+1; j < inst->nnodes; j++ )
{
sprintf(cname[0], "SEC on node pair (%d,%d)", i+1, j+1);
index[0] = xpos_compact(i,j,inst);
value[0] = 1.0;
index[1] = xpos_compact(j,i,inst);
value[1] = 1.0;
if ( CPXaddlazyconstraints(env, lp, 1, nnz, &rhs, &sense, &izero, index, value, cname) ) print_error("wrong CPXlazyconstraints on 2-node SECs");
}
}
...
*** SOME CPLEX'S PARAMETERS ***
// increased precision for big-M models
CPXsetdblparam(env, CPX_PARAM_EPINT, 0.0); // very important if big-M is present
CPXsetdblparam(env, CPX_PARAM_EPRHS, 1e-9);
CPXsetintparam(env, CPX_PARAM_MIPDISPLAY, 4);
if ( VERBOSE >= 60 ) CPXsetintparam(env, CPX_PARAM_SCRIND, CPX_ON); // Cplex output on screen
CPXsetintparam(env, CPX_PARAM_RANDOMSEED, 123456);
CPXsetdblparam(env, CPX_PARAM_TILIM, CPX_INFBOUND+0.0);
CPXsetintparam(env, CPX_PARAM_NODELIM, 0); // abort Cplex after the root node
CPXsetintparam(env, CPX_PARAM_INTSOLLIM, 1); // abort Cplex after the first incumbent update
CPXsetdblparam(env, CPX_PARAM_EPGAP, 1e-4); // abort Cplex when gap below 10%